function [OF, psi,rhoj,sj,s0j,etaEj,etaDj,vdj, piDprimej, dealer_ids_out, uB, uE, piB, piE, YB, YE, Ecost, uE2, psi_opt, rhoj_opt, W, Wxij, Wvdj, W2] = model_prediction_avg(q, Yt, params, side, SQ, CF, access,vdjSQ,idxj, t)

global Imax Dmax  xij_for_investor 

%sort Yt on id and then dealer 
Yt = sortrows(Yt, [1,16,17]);  %Note: this is done in the mainfile as well, should not do anything

if CF==0 && SQ==0
    estimate=1;
else
    estimate=0;
end

%index for dealers
dealer_ids  = Yt(idxj,16);
idxd        = find(ismember(Dmax,dealer_ids));

dealer_ids_out =  NaN*ones(Imax,1);
dealer_ids_out(idxd) =dealer_ids;

if CF>=8
q=q(idxd);
end 

%label params
sigma   = Yt(1,15);
munu    = params(1);
sigmanu = params(2);
c       = params(3);

%if CF==0
%   c=params(3);
%elseif CF==1
%    c=0;
%end

xij         = Yt(idxj,14);

if xij_for_investor==1
    xid     = 0;               %set to 0 if psi does not include xid
    
elseif  xij_for_investor==0
    xid     = xij;             %set to 0 if psi does not include xid
end 


%distribution functions
f = @(x) normpdf(x,munu,sigmanu);
F = @(x) normcdf(x,munu,sigmanu);
   

%initialize so that all outputes have the same number of dealers
psi        =  NaN*ones(Imax,1);
rhoj       =  NaN*ones(Imax,1);
psi_opt    =  NaN*ones(Imax,1);
rhoj_opt   =  NaN*ones(Imax,1);
s0j        =  NaN*ones(Imax,1);
sj         =  NaN*ones(Imax,1);
etaEj      =  NaN*ones(Imax,1);
etaDj      =  NaN*ones(Imax,1);
piDprimej  =  NaN*ones(Imax,1);
vdj        =  NaN*ones(Imax,1);
uB         =  NaN*ones(Imax,1);
uE         =  NaN*ones(Imax,1);
uE2        =  NaN*ones(Imax,1);
piB        =  NaN*ones(Imax,1);
piE        =  NaN*ones(Imax,1);
YB         =  NaN*ones(Imax,1);
YE         =  NaN*ones(Imax,1);
Ecost      =  NaN*ones(Imax,1);
mass       =  ones(Imax,1);

%use mean theta per dealer (to have a single theta per day)
theta= ones(Imax,1)*mean(Yt(idxj, 5));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Predictions for institutional
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%BUYERS  
if access==1 && side==1

        %psi
        psi(idxd,:) =  sigma*log(sum(exp(1/sigma*(xij+q))))-xid;
        psi_opt(idxd,:) =  max(xij+q);

        %market shares
        rhoj(idxd,:)     = mass(idxd).*(1-F(psi(idxd,:) +c- theta(idxd)));
        rhoj_opt(idxd,:) = mass(idxd).*(1-F(psi_opt(idxd,:) +c- theta(idxd)));

        %for non-benchmark dealer (deltajt=0)
        s0j(idxd,:)  = exp(1/sigma*(xij + q)) /(sum(exp(1/sigma*(xij + q))));
        %including bilateral
        sj(idxd,:)   = s0j(idxd,:).* sum(mass(idxd).*F(psi(idxd,:) +c - theta(idxd)));

        if CF==8 || CF==9 || CF==10 || CF==11
            s0j(idxd,:)=0;
            s0j(max(xij+q)==xij+q,:)=1;
            sj(idxd,:)= s0j(idxd,:).* sum(mass(idxd).*F(psi(idxd,:) +c - theta(idxd)));
        end 
        
        %elasticities 
        etaEj(idxd,:) = (1/sigma .* (1-s0j(idxd,:)) + sum(mass(idxd).*f(psi(idxd,:) +c  - theta(idxd)))./sum(mass(idxd).*F(psi(idxd,:) +c -theta(idxd))) .*s0j(idxd,:) ).*q;
        etaDj(idxd,:) = - mass(idxd).*f(psi(idxd,:) +c - theta(idxd))./(mass(idxd).*(1-F(psi(idxd,:) +c -theta(idxd)))) .*s0j(idxd,:).*q;
       
        if estimate==1
            %recover value of dealer
            a0(idxd,:)      = etaDj(idxd,:)./etaEj(idxd,:) .*rhoj(idxd,:)./sj(idxd,:);
            vdj(idxd,:) 	= (q .*(1+1./etaEj(idxd,:))+ (psi(idxd,:) +c- xij+xid).*a0(idxd,:))./(1+a0(idxd,:));
        elseif estimate==0
            %use the estimated value of the dealer
            vdj(idxd,:)     =vdjSQ(:,idxd)';
        end 

        %marginal profit 
        piDprimej(idxd,:) = (+xij-xid - psi(idxd,:) -c + vdj(idxd,:) ) .* etaDj(idxd,:) .* rhoj(idxd,:)./q;
        OF0 = sum((vdj(idxd,:) - q.*(1+1./etaEj(idxd,:).*(1-piDprimej(idxd,:)./sj(idxd,:)))).^2);

        if isnan(OF0)==0 
            OF = sum((vdj(idxd,:) - q.*(1+1./etaEj(idxd,:).*(1-piDprimej(idxd,:)./sj(idxd,:)))).^2);
        else 
            OF = 10^100;
        end 
        
        %EXPECTED PROFIT OF INVESTOR       
        %a) from bialteral trade
        uB(idxd,:) =rhoj(idxd,:).*xid;
        
        %b) from platform trade
        a       = psi(idxd,:) - theta(idxd) + c;
        Mills   = -f(a)./F(a);  
        ExpfunT = munu + sigmanu.*Mills ;
        uE(idxd,:)  = -ExpfunT + (1-rhoj(idxd,:)).*(-theta(idxd)+ (psi(idxd,:)+xid) );
        uE2(idxd,:) = -ExpfunT + (1-rhoj(idxd,:)).*(-theta(idxd) + sum((s0j(idxd,:).*(q+xij)))); %exclude sigma*epsilon
        Ecost(idxd,:) =  (1-rhoj(idxd,:)).*c;

        %Expected yield
        if F(a)~=1     %some trade bilateral
            YB(idxd,:) = rhoj(idxd,:).*(theta(idxd)-xij + xid) + (munu + sigmanu.*(f(a)./(1-F(a))));
        elseif F(a)==1 %all trade on the platform
            YB(idxd,:)=0;
        end    
        YE(idxd,:) = (1-rhoj(idxd,:)).* sum((s0j(idxd,:).*q)); 
        
        %EXPECTED PROFIT FOR THE DEALER
        %Profit from bilateral trades 
        piB(idxd,:) = (vdj(idxd,:) -theta(idxd) +xij - xid).*(1-F(a)) ...
                    - (munu + sigmanu*(f(a)./(1-F(a))));
                
        %Profit from platform trades
        piE(idxd,:) = sj(idxd,:).*(vdj(idxd,:)-q);
        
        %WELFARE
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)+xij-(theta(idxd)+munu)));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)));
        
        %Welfare when shrinking dealer's value towards theta
        distance = vdj-theta;
        vdj2     = theta+(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));
              
        
elseif access==0 && side==1 
    
        %use the estimated value of the dealer
        vdj(idxd,:)     =vdjSQ(:,idxd)';
    
        %Probability to trade OTC
        rhoj(idxd,:)= ones(size(idxd,2),1); 
        sj(idxd,:)  = zeros(size(idxd,2),1);  
        s0j(idxd,:)  = zeros(size(idxd,2),1);  
        
        %EXPECTED PROFIT OF INVESTOR    
        uB(idxd,:) =xid; 
        uE(idxd,:) =0.*xid; 
        Ecost(idxd,:)=0.*xid; 
        YE(idxd,:) =0.*xid; 
        YB(idxd,:) = theta(idxd) + munu -xij + xid;
 
        %EXPECTED PROFIT FOR THE DEALER
        %Profit from bilateral trades 
        piB(idxd,:) = vdj(idxd,:) -theta(idxd) -munu +xij - xid;
        %Profit from platform trades
        piE(idxd,:) = 0.*vdj(idxd,:);
        OF=NaN;
        
        %WELFARE
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)+xij-(theta(idxd)+munu)));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj(idxd,:)));
                      
        distance = vdj-theta;
        vdj2     = theta+(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));
       
    
%SELLERS    
elseif side==-1 && access==1
    
        %psi
        psi(idxd,:) =  -sigma*log(sum(exp(1/sigma*(xij-q))))+xid;
        psi_opt(idxd,:) = -max(-q+xij);
    
        %market shares
        rhoj(idxd,:) = F(psi(idxd,:) -c- theta(idxd));
        rhoj_opt(idxd,:) = F(psi_opt(idxd,:) -c- theta(idxd));

        %for non-benchmark dealer (deltajt=0)
        s0j(idxd,:)  = exp(+ 1/sigma*(xij - q)) /(sum(exp(+1/sigma*(xij - q))));
        %including bilateral
        sj(idxd,:)   = s0j(idxd,:).* sum((1-F(psi(idxd,:) -c - theta(idxd))));

        %elasticities 
        etaEj(idxd,:) = -(1/sigma .* (1-s0j(idxd,:)) + sum(f(psi(idxd,:) -c  - theta(idxd)))./sum((1-F(psi(idxd,:) -c -theta(idxd)))) .*s0j(idxd,:) ).*q;
        etaDj(idxd,:) = + f(psi(idxd,:) -c - theta(idxd))./(F(psi(idxd,:) -c -theta(idxd))) .*s0j(idxd,:).*q;
          
      if CF==0 && SQ==0
            %recover value of dealer
            a0(idxd,:)      = etaDj(idxd,:)./etaEj(idxd,:) .*rhoj(idxd,:)./sj(idxd,:);
            vdj(idxd,:) 	= (q .*(1+1./etaEj(idxd,:))+ (psi(idxd,:) -c + xij-xid).*a0(idxd,:))./(1+a0(idxd,:));
        else 
            %use the estimated value of the dealer
            vdj(idxd,:)     =vdjSQ(:,idxd)';
      end 

        %marginal profit 
        piDprimej(idxd,:) = (+xij-xid + psi(idxd,:) -c - vdj(idxd,:) ) .* etaDj(idxd,:) .* rhoj(idxd,:)./q;

        %OBJECTIVE FUNCTION 
        OF0 = sum((vdj(idxd,:) - q.*(1+1./etaEj(idxd,:).*(1+piDprimej(idxd,:)./sj(idxd,:)))).^2);
        if isnan(OF0)==0 
            OF = sum((vdj(idxd,:) - q.*(1+1./etaEj(idxd,:).*(1+piDprimej(idxd,:)./sj(idxd,:)))).^2);
        else 
            OF = 10^100;
        end 
   
        %EXPECTED PROFIT OF INVESTOR    
        %a) from bialteral trade
        uB(idxd,:) =rhoj(idxd,:).*xid;
        
        %b) from platform trade
        a       = psi(idxd,:) - theta(idxd) - c;
        Mills   = f(a)./(1-F(a));  
        ExpfunT = munu + sigmanu.*Mills ;
        uE(idxd,:)  = ExpfunT + (1-rhoj(idxd,:)).*(theta(idxd)- (psi(idxd,:)-xid)); %include sigma*epsilon
        uE2(idxd,:) = ExpfunT + (1-rhoj(idxd,:)).*(theta(idxd)- sum((s0j(idxd,:).*(q-xij)))); %exclude sigma*epsilon
        Ecost(idxd,:) =  (1-rhoj(idxd,:)).*c;

        %expected yield
        if F(a)~=0 %some trades billaterally
            YB(idxd,:) = rhoj(idxd,:).*(theta(idxd) +xij - xid) + (munu - sigmanu.*(f(a)./F(a)));
        elseif F(a)==0 %no one trades billaterally
            YB(idxd,:)=0;
        end 
        YE(idxd,:) = (1-rhoj(idxd,:)).* sum((s0j(idxd,:).*q)); 
         
        
        %EXPECTED PROFIT FOR THE DEALER 
        %Profit from bilateral trades 
        piB(idxd,:) = (theta(idxd) +xij - xid - vdj(idxd,:)).*F(a) ...
                    + (munu - sigmanu*(f(a)./F(a)));
                
        %Profit from platform trades
        piE(idxd,:) = sj(idxd,:).*(q-vdj(idxd,:));
              

        %WELFARE
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*((theta(idxd)+munu)-vdj(idxd,:)+xij));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(-vdj(idxd,:)));
              
        %Shrinking dealer's value towards theta
        distance = -vdj+theta;
        vdj2     = theta-(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));
       
    
elseif side==-1 && access==0  
    
        %use the estimated value of the dealer
        vdj(idxd,:)     =vdjSQ(:,idxd)';
        
        %Probability to trade OTC
        rhoj(idxd,:)= ones(size(idxd,2),1);  
        sj(idxd,:)  = zeros(size(idxd,2),1);  
        s0j(idxd,:)  = zeros(size(idxd,2),1);  
        
        uB(idxd,:) =xid; 
        uE(idxd,:) =0.*xid; 
        Ecost(idxd,:) =0.*xid; 
        YE(idxd,:) =0.*xid; 
        YB(idxd,:) = theta(idxd) + munu + xij-xid;
        
        %Profit from bilateral trades 
        piB(idxd,:) = theta(idxd) + munu +xij - xid - vdj(idxd,:);
        %Profit from platform trades
        piE(idxd,:) = 0.*vdj(idxd,:);
        OF=NaN;
        
        %Welfare
        W=sum((rhoj(idxd,:)+sj(idxd,:)).*((theta(idxd)+munu)-vdj(idxd,:)+xij));
        Wxij=sum((rhoj(idxd,:)+sj(idxd,:)).*(xij));
        Wvdj=sum((rhoj(idxd,:)+sj(idxd,:)).*(-vdj(idxd,:)));

        %Shrinking dealer's value towards theta
        distance = -vdj+theta;
        vdj2     = theta-(1-0.16)*distance;
        W2=sum((rhoj(idxd,:)+sj(idxd,:)).*(vdj2(idxd,:)+xij-(theta(idxd)+munu)));
        
end
